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Preheating after inflation involves large, time-dependent field inhomogeneities, which act as a 
classical source of gravitational radiation. The resulting spectrum might be probed by direct de- 
tection experiments if inflation occurs at a low enough energy scale. In this paper, we develop a 
theory and algorithm to calculate, analytically and numerically, the spectrum of energy density in 
gravitational waves produced from an inhomogeneous background of stochastic scalar fields in an 
expanding universe. We derive some generic analytical results for the emission of gravity waves by 
stochastic media of random fields, which can test the validity/accuracy of numerical calculations. 
We contrast our method with other numerical methods in the literature, and then we apply it to 
preheating after chaotic inflation. In this case, we are able to check analytically our numerical 
results, which differ significantly from previous works. We discuss how the gravity wave spectrum 
builds up with time and find that the amplitude and the frequency of its peak depend in a relatively 
simple way on the characteristic spatial scale amplified during preheating. We then estimate the 
peak frequency and amplitude of the spectrum produced in two models of preheating after hybrid 
inflation, which for some parameters may be relevant for gravity wave interferometric experiments. 



I. INTRODUCTION 



An expanding Universe with matter, radiation and dark energy is practically transparent for gravitational waves 
which, once produced, propagate freely to us. Cosmological backgrounds of gravity waves (see [l[ for early work, and 
0, 3, [3 for reviews) may thus carry unique and "clean" information about the universe at very early times / high 
energies. In particular, gravity waves play an important role in the context of inflationary cosmology. During inflation, 
tensor modes are produced from the amplification of initial quantum fluctuations into classical perturbations outside 
the Hubble radius, due to the accelerated expansion of the universe [f|. The resulting spectrum extends over a very 
large frequency range and its amplitude depends directly on the energy scale during inflation. These gravitational 
waves lead in particular to the B-mode polarization of the CMB anisotropy fluctuations. In models of inflation with 

large energy scales ~ (l0 15 GeV) (like chaotic inflation models), this imprint on the CMB anisotropies should be 
detectable by forthcoming CMB polarization experiments. On the other hand, if inflation occurs at lower energies (as 
in many hybrid inflation models), the amplitude of the resulting gravity waves would be too weak to be observed. 

However, there is another channel of gravity wave production related to inflation that is not associated with super- 
Hubble amplification of initial quantum fluctuations, but instead with the classical emission of sub-Hubble gravity 
waves from energy sources involving large, time-dependent inhomogeneities. Indeed, at the end of inflation, the 
inflaton decays and reheats the universe. In most models, the first stage of this process - preheating - is dominated by 
an explosive and non-perturbative production of highly inhomogeneous, non-thermal fluctuations of the inflaton and 
other bosonic fields coupled to it. In chaotic inflation models, the inflaton decays via parametric resonant particle 
creation Q , accompanied by violent dynamics of non-linear inhomogeneous structures of the (scalar) fields Q ■ In 
hybrid inflation models, the inflaton decays into inhomogeneous structures via tachyonic preheating Q, Q, [lol ]. 
The subsequent dynamics is characterized by turbulent interactions between Bose (scalar) waves, before the system 
ultimately settles into thermal equilibrium. The inhomogeneous decay of the inflaton, the subsequent turbulent stage 
and even the thermal state are inevitably accompanied by the emission of gravitational waves. The theory of gravity 
wave production from random media of inhomogeneous scalar field(s) after inflation has to address the question of the 
amplitude and typical frequencies of the resulting stochastic background. For low energy scale inflationary models, 
the frequencies of the gravity waves produced after inflation may well occur in the range which in principle can be 
detected by direct detection experiments (like LIGO/ VIRGO and BBO/DECIGO). This would provide us with a 
channel for testing inflation that can complement the CMB data, and with a unique observational window into the 
subsequent dynamics of the very early universe. 

There is no ready form of the theory of gravitational wave emission from random media that could be immediately 
applied for preheating. The formalism in Weinberg's textbook , based on the wave-zone approximation for localized 
sources in Minkowski space-time, is widely used. The energy in gravity waves is then expressed in terms of the double 
Fourier transform (in both space and time) of the stress-energy tensor. In principle, this formalism is applicable only 
for isolated sources, without expansion of the universe, and it does not allow one to follow the evolution with time of 
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gravity wave emission. There are several papers on the production of gravitational waves from a cosmological phase 
transition (see [l2j and references therein), due to the collision of bubbles and to the resultingdynamics. Gravity 
waves emitted from bubble collisions have been studied in [l3[ , [3] , [H[ , using the formalsim of [ll[ for the process of 
an instantaneous (i.e. fast compared to the expansion rate of the universe) phase transition. Calculations of gravity 
waves from hydrodynamical turbulence (l6j . |17j . departed from the formalism of [llj and solved instead for the 
wave equation of the Fourier tensor modes sourced by the transverse-traceless part of the stress-energy tensor. The 
considerations were restricted to the case of hydrodynamical incompressible fluids, and again for an instantaneous 
phase transition. 

The emission of gravitational waves from preheating was addressed in the original paper [l8j , and more recently in 



19(, [20( and [2l|. Refs. [18| and [19( performed numerical calculations of the gravity waves produced from preheating 



after chaotic inflation, using the formalism of For chaotic inflation models, the spectrum of gravity wave energy 
density has a bump at very high ferquencies, of order 10 s Hz, which is not likely to be observable. In 7|, two of us 
observed that the inflaton fragmentation in preheating after chaotic inflation occurs through the formation of non- 
linear bubble-like field inhomogeneities. It was conjectured in [3] that the main contribution to gravity waves emission 
comes from these bubble-like structures, both in the resonant and the tachyonic cases. Their relevance for gravity 
wave production in the case of tachyonic preheating was also pointed out in [10( . These bubbles have nothing to do 
with phase transitions. They originate as non-linear structures from the realization of the Gaussian random held of 
the vacuum fluctuations. If R is the charactersitc bubble size and H is the Hubble parameter when they form, then 
the fraction of energy density in gravity waves at the time of production is p gw /ptot ~ (RH) 2 at frequencies / ~ 1/R. 
Ref. [2(| studied numerically the production of gravity waves from parametric resonance at different energy scales, 
solving for the equations of the metric perturbations in Fourier space. Finally, Ref. plj introduced another numerical 
method to study gravity waves from tachyonic preheating after hybrid inflation. In particular, they displayed a case 
which, although involving very small coupling constants, they concluded could be detectable by BBO. 

In this paper we develop a theoretical framework for calculating systematically the emission of classical gravitational 
waves from random media of dynamical scalar fields in an expanding universe. The field inhomogeneities have very 
large (non- linear) amplitudes at scales smaller than the Hubble radius l/H but the media is homogeneous at larger 
scales. Our purpose is to extract observables like the present day energy density of gravity waves through statistical 
averages over the random scalar fields. We derive some general results for gravity waves produced by stochastic media 
of scalar fields, which are applicable for preheating, phase transitions and other instances of gravity wave emission by 
Bose fields. These general results will serve as tests for numerical methods suggested in earlier papers and as a check 
of our numerical results. 

We have used our formalism to develop a numerical algorithm for calculating gravity wave production. The field 
dynamics are calculated with lattice simulations, which are then used to find the evolution of the spectrum of gravity 
waves. These spectra are then converted into physical variables observed today. Our formalism gives significantly 
different results from others in the literature. Therefore we pay special attention to checking the results by analytical 
calculations. 

We illustrate our method with a model of preheating after chaotic inflation. This model provides a simple test-case 
for our formalism, allows easy comparisons with previous predictions, and allows us to verify our simulation results 
with detailed analytical calculations valid during the linear stage of preheating. 

The production of gravity waves from preheating after hybrid inflation is observationally more promising. We 
provide preliminary estimates for this case in this paper and will return to the issue in more detail in a subsequent 
publication. 

The paper is organized as follows. In section [II] wc develop the basic equations expressing gravity waves in terms of 
metric perturbations and relating these perturbations to a background of inhomogeneous scalar fields in an expanding 
universe. We also identify typical scalar field configurations that do and do not emit gravitational radiation. In section 
IIIII we compute the spectrum of energy density in gravity waves emitted by stochastic media, which includes taking 
spatial averages and ensemble averages over different realizations of the random fields. The spectra are then converted 
into present-day physical variables. In section HVi we outline analytical calculations of gravitational wave emission 
from different stages of preheating and thermalization, and we illustrate the formalism developed in the previous 
sections. We present our numerical method in section [Vj and we contrast it with previous numerical calculations of 
gravity waves from preheating. Section IVII contains our numerical and analytical results for the case of preheating 
after \<fi inflation. We discuss the influence of different stages of the scalar field dynamics on gravity wave production, 
and we express the amplitude and frequency of the spectrum's peak in terms of the model's parameters. We then 
perform a detailed analytical check of our numerical results, part of which is in an Appendix. In section IVII1 we 
estimate the amplitude and frequency of gravity waves produced in two different models of preheating after hybrid 
inflation. We conclude in section [Villi with a summary of our results and perspectives. 
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II. EMISSION OF GRAVITY WAVES BY INHOMOGENEOUS SCALAR FIELD SOURCES 

In this section we present the basic equations for classical gravitational waves emitted by a background made of 
inhomogeneous scalar fields. We also spell out a rather obvious but very useful theorem stating that the superposition 
of scalar fields waves with a "particle-like" dispersion relation does not emit gravity waves. 

A. Equations for Gravitational Waves with Sources in an Expanding Universe 

We consider several inhomogeneous scalar fields, denoted collectively by {<f> a , a = 1, 2, ...}, with energy-momentum 
tensor 

T^=d^ a d„<p a ~ gtiu (^gP° dpckad^a + V^j (1) 

where repeated indices a are summed. During preheating, the inflaton decays inhomogeneously and some of the fields 
coupling to it are significantly amplified. The fields are rather homogeneous at large scales (at the Hubble radius and 
beyond), but highly inhomogeneous inside the Hubble radius. These field inhomogeneities at small scales cannot be 
treated as small perturbations 1 . They participate in the evolution of the background scale factor through the average 
of the energy momentum tensor (T M „) in the Einstein equations. Here we consider the linear response of the metric 
perturbation Sg^ to the inhomogeneous part of T^. We will work at linear order in Sg^ because its coupling to 
T^ v is suppressed by the Planck mass Mpi, and the typical mass scales involved in the energy- momentum tensor are 
much lower than Mp\. Among the different components of the metric perturbation, gravitational waves are the only 
physical degrees of freedom which propagate and carry energy out of the source, see e.g. [22|. 

In a (spatially flat) Robertson- Walker background, gravitational waves may be represented by the transverse and 
traceless part of the spatial metric perturbation 

ds 2 = g^ dx^dx" = a 2 (r) [-dr 2 + (% + h tj ) dx'dx 1 ] (2) 

with 2 dihij — ha = 0. The perturbation hij corresponds to two independent tensor degrees of freedom and has the 
equation of motion 

h'l 3 +2^- ti tj - V 2 /iy = 16ttG a 2 ng T , (3) 

where a prime denotes a derivative with respect to conformal time t. 

Free gravitational waves obey the equation (|3|) without the source term. They can be quantized, and one can study 
the amplification of their vacuum quantum fluctuations in an expanding universe. An especially important case is 
exponential expansion of the universe during inflation, when tensor mode quantum fluctuations lead to a stochastic 
background of classical long-wavelength gravitational waves. 

We will consider a different, complementary situation, when quantum effects are negligible, and classical gravita- 
tional waves are generated by a non-zero source term in Eq. (J3]) . The source term n^ T is the transverse-traceless part 
(diHj^ — n^ T = 0) of the anisotropic stress Ily 

a 2 U i:j = Tij - (p) gij (4) 

where (p) is the background homogeneous pressure. Formally, the second term in the RHS of Eq. and the second 
term in the (i, j) components of the RHS of Eq. fl} involve the metric perturbation hij , through g^ = 8ij + hij . This 
gives a contribution in ^ (dk4> a dk4'a) hij (where hij is defined in ©) in the RHS of Eq. ([7]) below. We shall not 
include this term since it emerges only at second order in the gravitational coupling and is negligible at sub-Hubble 
scales. 

There are different ways to solve the wave equation ([3]). One can use Green functions in configuration space, but the 
transverse-traceless projection then involves inconvenient non-local operators. Another method [fl| . for the harmonic 



1 In fact, they correspond typically to density contrasts Sp/p ~ tens 

2 Here and in the following, Latin indices ... run over the 3 spatial coordinates, and repeated indices are summed. They are raised 
and lowered with the Kronecker symbol hij , so we don't distinguish between upper and lower indices. 
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gauge {d^hfiy = 0) in Minkowski spacetime, is based on the wave-zone approximation of the solution in configuration 
space, expressed in terms of the double Fourier transform (in both time and space) Tjj(k, ui) of the stress-energy 
tensor. This method was used in the context of preheating in Refs. [HI ]. [191 ], Formally, it is applicable only for 
isolated sources, without expansion of the universe, and requires the knowledge of the whole evolution of T^- with 
time. 

Here we develop another formalism, which is better suited to cases such as reheating with extended sources or 
continuous media in an expanding universe, and which allows us to follow the evolution of gravity waves with time. 
We will work in spatial Fourier space, with the convention 

r rfSu 

With the further redefinition 



Eq. © gives 



= a h i3 , (6) 



- h >>. (k ) + ( e — ) ha (k) = 16ttG a 3 n£ T (k) (7) 



where k 2 = k 2 is the square of the comoving wave-number. 

Given a symmetric tensor Ily , its transverse-traceless part is easily obtained in momentum space (but is non-local 
in configuration space) by the projection (see e.g. |23|) 



n£ T (k) = o Wm (k)n im (k) 



Pa (k) P jm (k) - - P ti (k) P lm (k) 



n im (k) (8) 



with 



P ij (k) = S ij -k i k j (9) 

where k = k/fc is the unit vector in the k direction. The operators Pjj are projectors on the subspace orthogonal to 
k, satisfying Pijki = and Pij Pji = Pu. From this it follows directly that ki n^ T = n? T = 0. 

In our case, the second term in the RHS of Eq. (j4|) and the second term in the components of the RHS 

of Eq. {TJ) are pure trace (gij = Sij, see the discussion below Eq. ((4])), and therefore they do not contribute to the 
transverse-traceless part. The relevant part of the energy-momentum tensor is then given by the product of the spatial 
derivatives of the fields, giving a convolution in Fourier space 

a 2 nJ T (k) = 7^ T (k) = 0^ lm (k) {dtfa d m <j> a } (k) = O ijllm $) J pi Pm Mp) <Mk - P) (10) 

where {di<f> a d m (f)a} (k) denotes the Fourier transform of di<f> a d m 4> a . In the last equality, we dropped a term in k m 
which vanishes when contracted with Ojj i j m (k). 

It should be clear that the process of gravity wave production that we consider here is very different from the one 
during inflation. During inflation, initial quantum fluctuations are amplified into super-Hubble stochastic classical 
perturbations, as a result of the accelerated expansion, that is through the term a" /a in Eq. ([7|)- In this case, the 
inhomogeneous part of the scalar field (inflaton) corresponds to a small perturbation (compared to the homogeneous 
part), and the source term n^ T vanishes at linear 3 order in 6<j)(x, r), see Eq. (jTU)) . This is to be contrasted to the 
case considered here, where the inhomogeneous part of the scalar fields cannot be treated as small perturbations, and 
n^.T act as a classical source. Note also that in Eq. ([7]), the term in a" /a vanishes when the equation of state is 
the one of radiation, w = 1/3. In general, this is not exactly satisfied during preheating. Indeed, when the inflaton 
potential is quadratic around its minimum, the average equation of state usually jumps during preheating from w = 
to an intermediate value that is close, but not exactly equal, to w = 1/3, see [27[. However, the dominant part 
of the gravity wave spectrum is produced well inside the Hubble radius at the time of production, so that a" /a is 



3 Gravitational waves can be produced classically from the second-order scalar perturbations generated during inflation ;24, 1251 . |26| . 
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negligible: a" /a ~ a 2 i/ 2 <C k 2 . In principle we can keep the term in a" /a and solve Eq. with a Green function 
more complicated than (fT2")) below. However, as we argued, this term can be dropped and Eq. ([7]) then reduces to 

ft£-(T,k) + fc 2 /iy(r,k) = 167rGa(r)^T( T;k) (n) 

where we have used Eq. (fTU|) . 

The solution with /lyC^i) = h!^(ji) — (assuming no gravity waves at the scale A: before the initial time Tj) is given 
by a simple Green function 

fty(T.k) = ^ f dr' sin[fc(r-r')] a(r') l£ T (r', k) (12) 



for k £ 0. Note that 

^ = 16ttG y dr' cos [fc (r - t')] a(r') 7^ t (t', k) . (13) 

If the source eventually becomes negligible after some time r = 17 (see below), the waves then freely propagate. The 
modes were sub- Hubble at the time of production and remain sub-Hubble until today. The corresponding solution of 
Eq. (jll[) without a source is simply 

hij(T, k) = Ay (k) sin [fc(r - r/)] + fly(k) cos [fc(r - r/)] for r > r/ . (14) 

Matching fry and TiJ • at r = 77 with (fT2l [13]) gives 

Ay (k) = 1^ jf ' dr' cos [fc ( r/ - t')] a(r>) T^(r', k) 

S«(k) = ^ ^dr' S m[fe(r/ -r')] a(r') ^ T (r',k) . (15) 

Before we develop our formalism further, we will discuss the emission of gravity waves by a special but interesting 
scalar field configuration. 

B. No-go Theorem: No Gravity Waves from Scalar Field Waves 



For the amplitude of the gravity wave, we obtained above 

^ 0«,l m (k) f dr' sin [k (r - r')] a(r') | ^ 



/iy(r,k) = — — O ijtlm (k) I dr' sin [k (r — r')] a(r') / /n _ V c, /9 Pm n (^, p) o (t', k - p) . (16) 



It is instructive to consider the very simple example of gravitational waves emitted by a medium made of a real free 
massive scalar field <fi obeying the Klein-Gordon equation in Minkowski spacetime (a(r) = 1 in the equations above). 
In this case, the Fourier transform of the field reads 

<Xp, r)e' ipx = 6(p) e ±4 ^ T+4 P x (17) 

with the dispersion relation 

u? v = p 2 + m 2 , (18) 

where m is the mass of <p. Decomposing the circular functions in Eq. (1 16|) or ()15|) into exponentials, and reordering 
the integrals over dr and d 3 p we then have 



Mr,k) oce ±ltT O y , im (k) J d 3 p PlPm b(p)b(k-p) J dT ' e i(±^±^| k - P i±fc)r' . ( 19 ) 

Eq. (|19[) corresponds to trilinear interactions between two </>-particles and one graviton, with different signs in the 
phase inside the time integral corresponding to different channels of interaction. In the limit of large time r compared 
to the frequencies of the particles, the time integrals reduce to Dirac delta functions enforcing energy conservation, 
such as Wp + wik-pi = k. Such trilinear interactions can be sketched by the diagram shown in the Figure[T] Momentum 



FIG. 1: Would be emission of a graviton hij with momentum k from the annihilation of two scalar waves 4>(p) an d <^(k — P) 
with momenta p and k — p. Helicity 2 of the emitted graviton cannot match the helicity zero of the incoming scalar waves. 

conservation is encoded in the convolution 4>{t\ p) </>(r', k — p). Energy and momentum conservations may be satisfied 
only for massless ^-particles, and only with k parallel to p. However, when k is parallel to p, the projector operator 
brings the gravity wave amplitude to zero, Oij,; m (k)p/p m = in Eq. (|19|). The reason for this is the conservation 
of helicity, which forbids interactions between free scalar waves and a graviton at linear order in the gravitational 
coupling. Interactions involving several gravitons are possible but further suppressed by the Planck mass. If instead 
of scalar waves the source of gravitational radiation is a superposition of vector field waves (photons or gauge bosons), 
gravitons can be emitted already at first order in the gravitational coupling, since diagrams like in Figure [1] but 
with vector fields carrying helicity 1 respect helicity conservation. For example, a thermal bath of photons emits 
gravitational waves [28[, while a thermal bath of massless scalars does not (see sub-section IIV Dl below) . 

Let us try to understand how far this no-go result can be extended, and which configurations of the scalar fields 
can lead to gravity wave emission. Suppose that instead of free scalar waves we deal with interacting scalar waves 
(which is typical for preheating). Then, instead of (fl8| . the dispersion relation of the field (f> will involve the other 
fields interacting with it. For instance, for a quartic interaction g 2 <jr tp 2 between </> and another scalar field ip, we 
have 

uj 2 = p 2 + m 2 + g 2 i\? . (20) 

If the frequencies uj p vary adiabatically with time, the time evolution of 4> can be described by the WKB approximation, 
4>(p,t) oc dr " p in Eqs. (|16p . (fT!?)) . Then trilinear interactions again satisfy energy conservation (technically, this 
results in this case from the stationary phase approximation), and the no-go result is intact. Therefore, we can 
formulate the no-go theorem: 

Scalar field configurations which can be represented as the superposition of waves with wave-like dispersion relations 
and adiabatically varying frequencies do not emit gravity waves at first order in the gravitational coupling. 

This result covers several interesting cases. For instance, during preheating after chaotic inflation, scalar field 
waves are produced only during short intervals of time and vary adiabatically between those instances. Therefore, we 
expect no gravity wave emission during the intermediate adiabatic regimes. We will treat this case in detail below, 
in sub-section IVI CI Between preheating and thermalisation the scalar fields enter a stage of Kolmogorov turbulence, 
characterized by the irreversible dynamics of weakly interacting scalar waves. We do not expect gravity waves to be 
emitted from this stage neither. Therefore, in the context of preheating, the time r = r/ introduced in the previous 
sub-section in Eq. (fTi]) . should be formally associated with the onset of the regime of weak turbulence. Another 
interesting case is a thermal bath of interacting scalar fields, which also does not emit gravity waves. We will discuss 
these cases in more details in Section HVl 

Now, let us try to understand which configurations of scalar fields can lead to the emission of gravitational waves? 
The no-go theorem is violated when the dispersion relation is not wave-like, and when the frequency does not vary 
adiabatically. We already mentioned one example of non-adiabatic change of the frequency of the field which is 
parametrically produced during preheating. In this case gravity waves are emitted even by massive scalars. Another 
example is the bubble- like configuration described by non-linear scalar fields where the wave-like dispersion is violated, 
as for bubbles formed from a first order phase transition. The collisions of scalar field bubbles emit gravity waves. 
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III. GW ENERGY DENSITY SPECTRUM 



In this section, we derive the spectrum of energy density in gravitational waves emitted by the random media of 
the scalar fields. This involves taking the spatial average of bilinear combinations of the (transverse-traceless) metric 
perturbation. This is what we will use in our numerical simulations. We develop this approach in sub-section IIII Al 
In sub-section IIII Bl we relate our formula to the Weinberg formula for -fjf 1 , which is often used in the literature. 
However, the simulations rely on a particular realization of the initial quantum fluctuations for the scalar fields. 
Analytical calculations can be advanced further by taking average over the ensemble of different realizations, and we 
calculate the gravity waves spectrum in this way in sub-section IIII CI Finally, in sub-section IIII Dl we present the 
rescalings needed to convert the spectrum into the present-day physical variables. 



A. GW Energy Density 



The energy density carried by a gravity wave cannot be localized in regions smaller than its wavelength, but it can 
be defined as an average over a volume V of several wavelengths' size (see e.g. (23l |) 



1 



Pgw 



32ttG 



(hij(t,ic) hij(t,x)) v 



(21) 



where a dot denotes a derivative with respect to cosmic time t = J a dr. In terms of conformal time and hij in ([6]), 
we have 



1 



hij hij — — h'ij + 2aH h t j h' tj + a 2 H 2 hij h lJf 



(22) 



For sub-Hubble wavelengths, k/a H, the second and third terms are negligible with respect to the first one. We 
therefore have 



1 



Pgw 



32nGa 4 



{KAt, x) h'iAr, x)) v 



1 



1 



32irGa 4 V 



d 3 k/4(r,k) h'*JT,k) 



(23) 



where * denotes the complex conjugate. In the second equality, we expanded each hij(r, x) into Fourier components 
hij(r, k), and then calculated the remaining spatial average as 



1 

V 



eTxe 



-i(k+k')x 



V 



<5 (3) (k + k') 



(24) 



where the (comoving) volume V has large dimensions compared to the (comoving) wavelengths A. The volume factor 
appears for dimensional reasons due to our use of a continuous Fourier transform ([5]), as opposed to a Fourier series. 
In the lattice simulations, V will correspond to the volume of the box in configuration space. The final results are 
independent of V. 

For hij (t, k) in Eq. (j2"3"|) , we use (fT4")) corresponding to the free waves propagating up to now after the emission 
process is completed. Suppose that today we are not interested in the resolution of the oscillation of hij(r, k) with 
time, so we average over a complete period of oscillation T = ^ 



(^(r,k)SS(r,k)) T = - £ (M + \Bi: 



(16ttG) 2 



E 

i,0 



T f 



dr' cos [k {r f - r')] a(r') T^ t (t', k) 



£ dr' sin [k (r f - t')} a(r') ^ T (r', k) ' j 



(25) 



where Y^i j = %ij X*j, and we used (fT5|) in the last equality. Expanding the cosine and sine above in factors of 

cos(fcT/) and sm(fcr/), these factors go out of the integrals and eventually give cos 2 (fcr/) + sin 2 (fcr/) = 1. Plugging 
the result into (|23|) . we get 



4ttG 1 



Pgw 



t,3 y 



dr' cos(fcr') a(T')^ T (r',k) 



dr' sin(fcr') a(r') T^ T (r', k) 



(26) 
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From l|26p , we construct the spectrum of gravity wave energy density per unit logarithmic frequency interval 

f dp gw \ = S k (T f ) 
\d\nk) T>Tf o*(r) 

where we have defined 

This is the main quantity that we will compute numerically, ft does not depend on the subsequent cosmological 
evolution, which dilutes the energy density in gravity waves and redshifts their frequencies. We will take care of the 
corresponding rescalings below in sub-section If 11 Dl The final spectrum today is then obtained from (|28p where t/ is 
the time when the source becomes negligible. However, we will also study how the gravity waves spectrum evolves 
with the time Tf before the source becomes negligible, in order to investigate the influence of the different stages of 
the scalar fields dynamics on the production of gravity waves. 

It is interesting to analyze Eq. (|28j) in the quadrupole approximation for large wavelengths (small k). In this 
approximation, one may drop the fc-dependence in the source termT^ T (r',k). We may then distinguish two different 
regimes where the spectrum of the energy density in gravity waves depends in a simple way on the frequency. For 
sufficiently small fc, the sine and cosine in the time integrals of Eq. (|28|) vary more slowly with time than the source 
and can be taken out of the integrals. In this case, the ^-dependence of S k comes only from the pre-factor in A; 3 , 
so the gravity wave spectrum varies as the cube of the frequency. On the other hand, for larger fc, there may be 
an intermediate regime where the quadrupole approximation is still valid but the sine and cosine vary more rapidly 
with time than the source. In this case, the time integrals of the factors in cos(fcr') and sin(fcr') in Eq. (|2"8)) are 
proportional to 1/fc. Overall, this gives S k oc k, so the gravity waves spectrum varies linearly with the frequency. 



dr cos (/cr) a(r')T^ T (r',k) 



dr' sm(fc-r') o(r') T^ T (r', k) 



(28) 



B. Relation to the Weinberg Formula for d ^ w 

The formalism developed above is designed for cases, such as reheating, where the source or continuous media 
extends over the whole expanding universe, and is active - with respect to gravity wave emission - for a limited period 
of time (until Tf). It may be directly extended to systems emitting gravity waves continuously - i.e. active all the 
time - by taking the limit Tf — > oo. In this case, we may relate our formalism to the one derived by Weinberg (l]| 
for isolated sources in the wave-zone limit (i.e. at distances large compared to the wavelengths and to the size of the 
localized source) in Minkowski spacetime, which is often used in the literature. 

Consider the double Fourier transform of the energy-momentum tensor in both space and time (in Minkowski 
spacetime) with the conventions of [ll[ 

Ta (k, w ) = J ^ e wr J d 3 x e- ikx T l3 (r, x) . (29) 

According to [ll[ , the total gravity wave energy per element of solid angle emitted by the localized source is given by 



dE gw 
dQ 



2G O iUm (k) J dk k 2 Tij (k, k) r; m (k, k) . (30) 



For future reference, we now outline how this formula can be derived in our formalism. First, we set a(r) = I and 
take the limit Tf — > oo in Eq. (|26p . Next, we develop the sine and cosine into exponentials, and express T? t (t, k) in 
terms of Tij(k, uj) with ([5]) and (|I0p . The resulting bilinear combination of T!y(k, u) involves a bilinear combination 
of the projection operators Qij t i m (k), which satisfy 

&ij,lm (k) Oij^rs (k) = C; m rs (k) . (31) 

Simplifying Eq. (|26[) for the gravity wave energy density accordingly, we recover Eq. (|30[) for the total energy in gravity 
waves per element of solid angle in k-space. 

The scale factors in our formula arise because, with the expansion of the universe, the energy density (|26[) dilutes 
as radiation, p gw oc a~ 4 , and the factor of a(r') in the integrals comes from the fact that r is the conformal time (and 
k the comoving wave-number). Once again, note that in (|28p we may follow the evolution of gravity wave emission 
with the time 77. 
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C. Ensemble Average 



So far we have considered the emission of gravity waves by inhomogeneous media, without using their stochastic 
character. We now consider the gravity wave energy density (f2"TT) obtained from the ensemble average over different 
realizations of the random scalar fields. Proceeding as in (f2"2")) and decomposing hij (r, x) into Fourier components, we 
take the ensemble average (...) of the resulting bilinear combination of the gravity wave amplitudes 



P S w = ^ZTTZZ I 7^. I (^,(k,r)^(k',r)) . (32) 



d 3 k r d 3 k' 

327rGa 4 J (2tt) 3 / 2 J (2tt) 3 / 2 

For hij(r, k), we take again the free waves (TT¥1) - (|15[) after the emission process is completed. This involves the 
calculation of 

(A;(k)A*.(k')> = <^¥- dr' fj dr" cos [k(r f - r')] cos [k'(r f - r")\ a{r')a{r") (T^(r', k) T^*(r", k')) . 

We are thus led to consider the unequal time correlators of the transverse-traceless part of the energy momentum 
tensor ([Toll 



<i;7(T',k)7;7V,k')) = 

O ij>lm (k)O ijirs (k>) fj0jl |^^p^ m p;p' s (^(p,r')0a(k-p,r')^(p',T'')^(k'-p',r'')) . (34) 

For an arbitrary random medium, this depends on the 4-point unequal time correlators of the scalar fields. We can 
proceed further by assuming that the <j) a {p) are described by a statistically homogeneous random Gaussian field. This 
description is often relevant in the context of preheating. The preheating process is conveniently divided into four 
different stages: linear preheating, non-linear preheating with sign ificance of backreaction, Kolmogorov turbulence of 
scalar waves, and thermalization. Numerical investigations [7, 29, 30] show that the scalar fields are Gaussian during 
the linear stage of preheating and non-Gaussian at the transient, "bubbly" stage and at the early turbulent stage, 
while gaussianity is recovered at the later stage of turbulence and at the thermal stage. The following analytical 
analysis can be implemented for fields with Gaussian statistics. The 4-point functions may then be developed in 
terms of 2-point functions (cf. the Wick theorem) as 

(^(p,r')0 Q (k-p,r')^(p',r")^(k'- P ',T")) = (^(p, r')^(p - k, r')) (&(-p', r") <j? b {k! - p', r")) 

+ (0 q ( P ,t')^( P ',t")) (<Mk-p,T')^(k'-p',T")) + (<Mpy)^(k'-pV)) (<Mk-p,T')C(pV)) (35) 

where we have used 4>* a (p) = 4> a {—p) for real scalar fields. The unequal time correlators of the scalar fields are 
expressed as 

(^(p,t')^(p',t")> = F ab (p,r',r") S(p-p') (36) 

where F ab depends only on p = |p|, by statistical homogeneity and isotropy of the scalar fields. The first term in 
the RHS of (|33|) is proportional to S(k) S(k') and will not contribute (to the connected part of the 4-point function). 
Using with in §M$, we then get 

(T^(r',k)T^(r",k')) = 

5(k-k')O ijM (k)O ij<rs (k) J ^3 [pip mPr p s +pip m (k r - Pr ) (k s - Ps )} F ab (p,T' ,t") F ab (\k - p\, T ' ,t") (37) 

where we have performed the integral over p' and enforced k = k'. In this expression, the terms proportional to k r 
and k s vanish when contracted with Oy irs (k). Using Eq. (|3I[) . we then have 



I r 



2 



.,4 



Ol m ,rs{^)PlPrnPrPs = ^ p 2 - (k.p) 2 = ^- SU1 4 (k, p) (38) 



2 



where (k, p) denotes the angle between k and p. Simplifying (|3"7| accordingly, and inserting the result into (f3"3")> . we 

get 

(A l3 (k) A% (k')) = <5(k - k') (16 ^ 2 G) £ dr 1 £ dr" cos [k(r f - r')] cos [k( Tf - /')] a(r') a(r") 

/ 7^P 4sin4 ( k 'P)^fe T '^'')^(|k-p|,r',r'') (39) 
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Averaging over a complete period of oscillation of the free waves as in (l25t , we then have 

((^■(k, r) 4*(k', r))) r = y «^(k) A^k')) + (B«(k) B*.(k')» (40) 

where (Bjj(k) S*-(k')) is given by (|39[) with the two cosines replaced by sines. These combine in ([40]) to give 
cos [k(r' — t")]. Finally, collecting all the factors in (f3"2")l . and defining the gravity wave energy density spectrum 
as in l[S7j). we find 

S*(r/) = \ G k 3 J j^L p 4 sin 4 (k, p) f dr' J*' dr" cos [k{r' - r")} a(r') a(r") F ah {p, r', r") F ab (\k - p|, r', r") 

(41) 

where we have performed the integral over k' in (|32[) . as well as the one over the direction k. This last step comes 
from the fact that the result of the integral over d 3 p in (|4Tj) is independent of the direction k, which is to say the 
gravity wave spectrum is isotropic. 

Eq. (|4"Tj) is the main result of this section. We will use it in Section UVl and in sub-section I VI Cl to perform analytical 
calculations of the gravity wave energy density spectrum. 



D. Spectrum Today 

Eventually we will be interested in the abundance of gravity wave energy density today 

Pgw \ _ f df u2 , 



= / ^fi gw (/), (42) 

Pc J o J J 



and its spectrum per logarithmic frequency interval 



where / is the frequency and p c = 3Hq / (8ttG) is the critical energy density today. To convert the gravity wave 
spectrum considered in the previous sections into physical variables today, we need to consider the evolution of the 
scale factor from preheating up to now, which depends on the evolution of the equation of state. In general, when 
the inflaton potential is quadratic around its minimum, the equation of state jumps from w — to an intermediate 
value close to w = 1/3 during preheating (27j . (3l| . We assume that the universe then evolves continuously towards 
radiation domination (e.g. without an intermediate matter dominated stage ). 

Let us denote by ti the end of inflation, by tj a moment after the jump of the equation of state (to be specified 
below), by > tj the moment when thermal equilibrium is established, and by to the present time. From t„ to to we 
assume the expansion of the universe obeys entropy conservation, and from tj to t* we assume that it evolves with 
a mean equation of state w (which is thus close to 1/3). The scale factor today compared to the one at the end of 
inflation may then be expressed as 

/ \ 1-§(1+M>) / \ -1/12 

°±-^(«l\ (44) 



a o d, /9- /4 \ a *J \9o 



where pj is the total energy density at t — tj, p rac io is the energy density of radiation today, and g is the number 
of effectively massless degress of freedom (neglecting the difference between g and gs), see [32|]. From Eq. (|4*4")h we 
deduce the physical wave- number today, ko = k/ao, in terms of the the comoving wave- number k used in the previous 
sections. The corresponding physical frequency today is given by 

kn k /n V^ 1 ^ 



4 An intermediate stage of matter domination could be due, for instance, to a massive relic dominating the expansion before it decays. If 
such a stage starts when the scale factor is a(t) = Oi, before the universe becomes radiation dominated at a(t) = 0,2, then compared to 
the relations given below, the present-day gravity wave frequency / would be smaller by a factor of (ai/c^) 1 / 4 , and the amplitude of 
the spectrum f2 gw h 2 would be smaller by a factor of (ai/cta)' 
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The energy density spectrum today is given by ([27]) with a = do- Using (|44[) . this gives 

, T . „ , n n 1 —371! , - -1 /a 



where £l la d h 2 = /i 2 p ra do/Pc = 4.3 x 10 -5 is the abundance of radiation today, and we will take g*/go = 100. 

In the above relations, dj and pj may be derived directly form the simulations for a precise calculation. Everything 
is then known, except for the factor a^/a*, which depends on the unknown reheating temperature. Note that this 
factor reduces to one for w = 1/3. Therefore, by chosing tj such that Wj is already close to 1/3, this factor gives only 
a small contribution (which may be bounded by imposing some constraints on the reheat temperature T*). 

In the case of a A </> 4 model of inflation, the average equation of state reaches w — 1/3 already during the simulation, 
in fact soon after the end of inflation, so the factor in a^/a* is just absent in this case. The above relations then 
reduce to 

For A (f> 4 model: / = k 4 x 10 10 Hz and ft gw h 2 = 9.3 x 10~ 6 ^fll , (47) 

where tj will be chosen at the end of the simulations. Since in this model the average equation of state is the one of 
radiation very quickly after the end of inflation, we have aj pj ~ a\ pi, so that, normalizing the scale factor to one at 
the end of inflation, aljpj is essentially given by the energy density at that time, a^pj ~ A 0^/4 ~ 3.4 x 10~ 17 Mpj (for 
A HI ; >. 



IV. ANALYTICAL CALCULATIONS OF GW EMISSION 



In this section we outline analytical calculations of the gravity wave spectrum emitted from preheating, and we 
illustrate the general formalism developed in Sections [Til an d IIIIl We consider four different stages of reheating: the 
linear preheating stage, the non-linear "bubbly" stage, scalar wave turbulence, and the stage of thermal equilibrium. 



A. Linear Preheating Stage 



Consider the first, linear stage of preheating when only one field, call it Xi nas been amplified. The quantum field 
X is described by the field operator 



d 3 k 

(2tt) 3 / 2 



(48) 



where and dk are creation and annihilation operators with usual commutation relations. The field x obeys 
Gaussian statistics, so that the formalism of Section Till CI should be applicable here. In Eq. (|36|) . instead of the 
Fourier amplitude x(p,r), we have to use the operator 5 



The field correlator is 



<0|x(p, r') X + (p', r")|0) = Xp (r') X^") S(p p') 



(49) 



(50) 



which gives F xx (p, t', t") = Xp( T ') Xp( T ") m ([3^| . Inserting this into Eq. (|4T|) . and expanding the cosine, the integrals 
over t' and r" are then complex conjugates of each other 



Skirt) 



Gk : 



^ 3 P 4 • iff -\ 

p sin (k, p) 



(27T) 



dr cos (fcr) a(r) x P (r) X|k- P | (r) 



dr sin (k r) a(r) Xp(r) X|k- P | ( 7 



(51) 



5 Note that the modes \ p in II49I I have dimensions mass 1 / 2 , while x(P> r ) has dimension mass 2 . 
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This equation allows us to perform analytical calculations of the gravity wave spectrum emitted in different models 
during the linear stage of preheating, where the time evolution of the eigenmodes Xp( T ) may be derived analytically. 
For example, in the case of broad parametric resonance with a quartic interaction g 2 (j) 2 x 2 > the dispersion relation for 
the Xp waves is 

uj k = {k/af + g 2 4> 2 (t) , (52) 

where <p(t) denotes the background inflaton oscillations at the end of chaotic inflation. For tfi away from its zeros, the 
frequency cu^ is changing adiabatically with time, and the Xp( r ) are described by the WKB solution. According to the 
theorem of the sub-section III B| no gravity waves are emitted during these stages. However, when <j>(t) crosses zero, 
u>k changes non-adiabatically and gravity waves are emitted. We will calculate in detail the gravity wave spectrum 
emitted from broad parametric resonance below in sub-section I VI Cl for the model (I66|) . As we will see, our numerical 
results are reproduced very accurately from the analytical theory. 

B. Non-linear Bubbly Stage 

At the end of preheating, there is a violent, highly non-linear and non-perturbative stage, where the inhomogeneous 
fields have very large occupation numbers and strongly interact. The fields then become strongly non-Gaussian and 
analytical forms for their amplitudes are not available. Therefore, the analytical method of the previous section is 
not very useful for this stage. Yet, we expect this stage to give a dominant contribution to the gravity wave emission. 
It will thus be very useful to have an analytical estimate of the effects. 

For preheating after chaotic inflation, visualization of the field dynamics in configuration space 1 reveals that, 
during the linear stage of preheating, scalar fields interacting with the inflaton are produced as standing random 
gaussian fields with exponentially increasing amplitude, then non-linear rescattering generates standing random non- 
gaussian inflaton inhomogeneities with very fast growing amplitude. The peaks of the inflaton inhomogeneities 
coincide with the peaks of the scalar fields produced by parametric resonance. When the inflaton peaks reach their 
maxima, they stop growing and begin to expand. The subsequent dynamics is characterized by the expansion and 
superposition of the scalar waves originated from the peaks. Multiple wave superposition results in phase mixing and 
sets up the turbulent waves dynamics. Thus, the short intermediate stage is characterized by the formation, expansion 
and collision of bubble-like fields inhomogeneities associated with the peaks of the Gaussian field. This process is 
qualitatively similar to the bubble-like inflaton fragmentation occuring in tachyonic preheating after hybrid inflation 
1,1,0. 

The typical (physical) size i?* of the bubble- like fields inhomogeneities when they fragment depends on the character- 
istic (comoving) momentum amplified by preheating, ~ a/fc*. It was conjectured in [7[ that this fragmentation 
gives the dominant contribution to gravity wave emission, with a total energy density in gravity waves of order 

fPsA ^a(R*H)l (53) 

VPtot/p 

at the time of production t p , where H is the hubble radius at that time (i.e. during the "bubbly" stage). We will 
estimate the coefficient a from numerical simulations in sub-section IVI Bl We will see that, for the model (|66D . 
Eq. (|53p reproduces well the numerical calculations of the gravity wave emission, for a ~ 0.15. Eq. (|53|) seems rather 
general. It already appeared in the context of gravity wave emission from the collision of bubbles formed from first 
order phase transitions (in the thin- wall approximation) [33| , [l3j . For comparison, in the case of the collision of two 
(relativistic) vacuum bubbles, Ref. [l3[ finds a = 1.3 x 10 -3 . 

C. Scalar Wave Turbulence 

After preheating, the scalar fields enter a stage of turbulent interactions. The well-developed turbulent stage 
corresponds to media of the scalar waves, weakly interacting and having random phases. They have dispersion 
relations similar to (|18[) and are described by eigenfunctions varying as e~ lu>kt with time. Thus, according to the no- 
go theorem of sub-section III Bl no gravity wave emission is expected from the scalar waves' Kolmogorov turbulence. 
There are, however, some subtleties. In fact, numerical simulations of the scalar field dynamics show that developed 
Kolmogorov turbulence is not established immediately after preheating ends. One of the indicators of this is that 
the fields stay non-Gaussian for a while after the end of preheating [3, H(|. Also, for models of chaotic inflation, 
the residual inflaton condensate <fi is still significant after the end of preheating [27|, Hi]]. Therefore, the wave-like 
dispersion relation (|18[) is not established immediately at this stage. We shall expect some residual gravity wave 
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production after the "bubbly stage" at the end of preheating. This is hard to treat analytically but can be seen in 
the numerical calculations of Section IVI1 



D. Thermal Bath 



Consider a scalar field <p in thermal equilibrium with a thermal bath at a temperature T . As long as the temperature 
is higher than the Hubble parameter H, the expansion of the universe is not relevant. The scalar field is described by 
the field operator decomposed into oscillators as in (|48p. which again obey Gaussian statistics. Here, instead of the 
Fourier transform of sub-section IIII Cl we have to use the operator 

<p(Kr) -» — =a k + -=ai k . (54) 
V ZuJk v^fc 

where uj\ = fc 2 + m 2 and m is the mass of ip. Then we can directly use Eq. (|4"l"j) , where we have to calculate the 
unequal time field correlator (f3"6"| . In order to do this, we take an average for the thermal bath using the density 
matrix Z~ x e~ H / T 

^(k,r)^(k',r')) =Z- 1 Tr( e -T^(k,T)^(k',T')) , (55) 

where the partition function Z = Hk(l — e~^) and the Hamiltonian for oscillators H = X^k^ka^ak- Next, we sub- 
stitute in (|55p the operators in the form (|54[) . Dropping the negligible contribution from the zero vacuum fluctuations, 
the result for the thermal, non-vanishing part is 

Z-'Tr (e-%a+au) = V(n p | e -^a+a k ,|n p ) = . (56) 

\ / g T — 1 

Then, for a scalar field in equilibrium with a thermal bath at temperature T, the unequal time field correlator is 

(y(p,r)yVy)) - C ° S[ <i {T /T T l *(p~pO - ( 57 ) 
cj p [e u pl 1 - 1} 

so that 

w p (e^p/- 1 — 1J 

in Eq. (|36j). 

The field correlator is a harmonic function of time. Therefore after inserting it in Eq. (I4ip . we again reproduce the 
conditions for the no-go theorem of sub-section III Bl No gravity waves are emitted from the thermal bath of scalar 
fields. However, the bath of gauge fields and photons emits gravity waves (28|. 



V. NUMERICAL METHODS FOR GW CALCULATION 



In this Section we describe our numerical algorithm, which is based on the formalism of Section IIII1 Then we 
compare our method with previous numerical calculations of gravity waves produced from preheating. 



A. Numerical Algorithm 

Our method is based on numerical computation of the quantity Sk{rf) defined in (|28|) . We use the program 
LATTICEEASY [35j to calculate the evolution of the scalar fields and the evolution of the scale factor. A typical 
realization of the scalar field profile in the most interesting "bubbly" stage (during preheating) is shown in the 
Figure [2] The calculation of the scalar field evolution is performed in configuration space and is described in the 
LATTICEEASY documentation. We added a function to periodically calculate the integrands in equation (|28l) . 

To find the integrand at each integration time step we calculated the transverse, traceless components of the scalar 
fields' energy momentum tensor T^ 1 (r' ,k), given by equation (|10[) . For each field we calculated six arrays with the 
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mixed derivatives {di4>dj<fi} (x) in configuration space, and then Fourier Transformed them to obtain the six arrays 
{9i^^}(k). 

We found the projection operator Oij.im analytically in each direction we used, using Eqs. ([8} and Q. For example, 
from Eq. (0) we can see that in the k x direction the only nonvanishing components of are P22 — P33 = 1- From 
Eq. ([8]) you can thus show that the only nonvanishing components atOijj, m along this axis are ©22,22 = 033,33 = (1/2), 
022,33 = 033,22 = —(1/2), 023,23 = 032,32 = 1- Finally, we can plug this into Eq. |T]) to find the only nonvanishing 
components of T? t (t' , k) along the k x axis: 

Tl 2 T = -Tj 3 T = i[{ft0fc#(k)-{fcW}(k)] (59) 
T^ = T^ = {&<^}(k) 
In a similar way one can derive expressions for the components of T^ T (r',k) along any direction in k space. 




FIG. 2: Two-dimensional slice through a three- 
dimensional realization of the scalar field \ from a numer- 
ical simulation of preheating in the \(f> 4 + g 2 (j> 2 X 2 model 
of the next section. The horizontal axes correspond to 
spatial coordinates and the vertical axis corresponds to 
the field's values. 



FIG. 3: Spectrum of energy density in gravity waves cal- 
culated along nine different directions in k-space. The 
parameters for this run were the same as for Fig. [4] see 
sub-section I VI Bl for details. 



We calculated the values of Sk along nine different directions in Fourier space. The first three were the three axes 
k x , k y , and k z and the remaining six were directions at 45° in between two of the axes: xy, xz, yz, {—x)y, (-x)z, 
and (— y)z. The results for the spectrum £l gw h 2 calculated along these nine different directions is shown in Figure [3] 
This figure shows that Sk is nearly identical in different k-directions (i.e. the gravity wave spectrum is isotropic), 
but it also shows that there is statistical noise in each spectrum plot. Therefore, we averaged over the six diagonal 
directions, which significantly reduces this noise. Unless otherwise noted all spectra shown in this paper represent 
averages over the six diagonal directions in k space. (We did not include the axes because in a cubic lattice the k 
values along the axes are different from the ones along diagonals, making it harder to average the results together.) 
In Figures 2] and [TU] an average over the three axial directions is plotted as well as an average of the six diagonal 
directions. 

Finally, we used Eq. (|4T|) to calculate / and fl gw h 2 today. We confirmed that the combination afipj was constant 
from a time shortly after inflation to the end of the simulation, as it should be for radiation domination, so we used 
the values from the end of the simulation. 



B. Comparison with Other Methods 



In the literature, there are different interesting methods of numerical calculation of gravity waves emission from 
preheating. In this sub-section, we contrast our formalism to previous numerical methods. 

Gravitational waves from preheating after chaotic inflation were investigated in [l8|], [l!|, on the basis of the 
Weinberg formula (|30[) in Minkowski spacetime. The expansion of the universe was then taken into account by 
dividing the conformal time into steps At, calculating the gravity wave energy density resulting from each step, and 
summing up each contribution diluted by the value of the scale factor at the middle of the step. In addition to 
the crude account of the evolution of the scale factor, one may show from sub-section IIII Bl that such a treatment 
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effectively replaces the terms 



in Eq. (J2SJ) by 



L 



dr' cos(fcT') a(r')^ T (r',k) 



(60) 



E 



,+At 



dr' cos(fcr') o(r')T^ T (r',k) 



(61) 



where a — 1,2,3, ... counts the steps. Indeed, what is summed up in Refs. [l8| , [l9j is the energy density from each 
partial step, instead of the energy density from the whole evolution of gravity waves with time. To compare this with 
our method, we can divide the time integral in (I60|) into similar steps. Obviously, the result is different from (|61l) . 
which misses all the cross terms given by the products of integrals over different steps. In other words, the approach 
based on (|6ip does not take into account the precise propagation history of the gravitational waves in the medium. 
Nevertheless, this may give a good approximation during preheating itself, when the scalar fields are exponentially 
amplified and their past values give negligible contributions. However, the past evolution should be taken into account 
at the end of preheating, when the amplitude of the scalar fields stabilize around their maximum value. As we will see 
in the next section, our numerical results are relatively similar to the ones obtained in the earlier paper [l8| (where 
a smaller frequency range is considered), but they differ significantly from the results of [l9j |. for both the shape and 
the amplitude of the gravity wave spectrum (see sub-section IVIB|) . 

Another method to calculate the production of gravitational waves from preheating was used recently in (2lj . in 
the context of hybrid inflation models, without expansion of the universe. Ref. [2l| considers the traceless part hij 
(ha — 0) of the spatial components of the linear metric perturbation in the harmonic (de Donder) gauge (d^h^ = 0), 
satisfying 



16ttG 



(62) 



and solve for this equation in configuration space. Gravity waves correspond to the two independent degrees of freedom 
of the transverse and traceless part of the spatial components of the metric perturbations, i.e. the tensor part, that 
we denote hy. However, in addition to the tensor part hij, the five independent degrees of freedom of the traceless 
hij in harmonic gauge involve an extra scalar part (with one degree of freedom) and an extra (divergence-less) vector 
part (with two degrees of freedom). In configurations space, one has to solve for these extra components of the metric 
perturbation as well, because the transverse-traceless projection ([5]) is then non-local. Only the transverse-traceless 
part contributes to the energy density in gravity waves, and therefore it has to be extracted efficiently 6 . Ref. [2l[ 
proposes to use the formula 



2 - - 

(p gw oc) (hij (x) hij (x)) = - (hij (x) h tj (x)) 



(63) 



to relate the energy density in gravity waves to the spatial average of all the components of the traceless hij . 

We can justify Eq. (|63[) in the special case of gravity waves emitted far away from an isolated source, in the 
quadrupole approximation (i.e. for wavelengths large compared to the size of the localized source). Indeed, Eq. (j6"2"|) 
can be solved in Fourier space with the Green function similar to the solutions 1]12[) . (flU)) (here in Minkowski spacetime) 



hij(t,\s) 



dt cos[fc(i-f')] Mij(t',k) , 



(64) 



where My (t, k) is the Fourier transform of the RHS of Eq. 



The transverse-traceless ha is then obtained from 



the traceless hij by the projection ([5]) in Fourier space, which satisfies (pTIj) . The spatial average of the bilinear 
combination of the tensor modes involved in the gravity wave energy density is then given by 

d 3 k 



(hij(t,x)hij(t,x)) = 



(2tt) 3 /' 



Or 

m,rs (k) him, (t,k) h rs {t,k) 



(65) 



In Weinberg's formalism in sub-section IHI Bl which is also derived in the harmonic gauge, the separation of the gravity waves from the 
other modes is achieved by the wave zone approximation far away from the source and by imposing the dispersion relation ui^ = |k| for 
the solution of the wave equation. 
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In the quadrupole approximation, the Fourier transform of the source, My (i,k), does not depend on k, so that 
hij(t, k) in (fM]) does not depend on the direction k. In this case, the integrand in the RHS of Eq. (|6"5"]) depends on 
the direction k only through C; m!rs (k), giving a 2/5 factor after integration over the solid angle in k-space. This then 
leads to Eq. (f6"3"]) . However, it is not clear how this equation is applicable outside the quadrupole approximation and 
for extended sources. 

Indeed we can give an example of scalar fields configurations relevant for preheating where the method of [2l| is not 
applicable. As in sub-section III B( consider media of superposed scalar field waves. There is a gravitational response 
to its inhomogeneities and dynamics in the form of scalar modes, but the tensor modes arc not emitted in this case, as 
we derived in sub-section III Bl Indeed for this case, in (j6"3"]l we have (hij(x) hij (x)) = 0, but in general (except in the 

quadruple approximation) (fty(x) /ijj-(x)} ^ 0. This may be checked explicitely with the solution (|64[) by proceeding 
as we did in sections III Bl and Mil The argument of sub-section III Bl does not apply for hij (in other words, trilinear 
interactions between hij and two scalar particles do not violate helicity conservation), because it involves a scalar 
degree of freedom. The method of [2l[ is thus not generic, since it would lead to the incorrect conclusion that there 
is a non- vanishing emission of gravity waves from the media of scalar waves. 



VI. APPLICATION TO PREHEATING IN THE A (j> A MODEL 

In this section, we illustrate our results for a simple model of preheating after A tfi 4 chaotic inflation. We first review 
some useful features of the scalar field dynamics. Then we present our numerical results for the spectra of energy 
density in gravitational waves. Finally, we give an analytical treatment of the gravity wave production during the 
linear stage of preheating, which allows us to check our numerical results in detail. 



A. The Model 



The model we consider corresponds to two scalar fields minimally coupled to gravity with the potential 

V=- A ^ + 9 ^ 2 X 2 (66) 

where A — 1CP 14 (from normalisation of CMB anisotropics) . The theory of preheating for this model was developed 
in [36| . Because of conformal invariance, the expansion of the universe may be scaled out of the equations of motion. 
It is convenient to redefine the fields and the time coordinate as 

X = a\ , ^ — acj) , x — VXcpoT (67) 

where 4>o ~ 0.342 Mpi is the initial amplitude of the inflaton condensate at the end of inflation, and the scale factor a 
is normalised to 1 at that time. During the linear stage of preheating (when backreaction is negligible), the temporal 
part of the modes Xk (x) obey the oscillator equation 7 

' /_ A k + uj 2 k (x) X k = with u\ (x)=K 2 + qp (x) (68) 



dx 2 

where / denotes the inflaton zero- mode rescaled by its initial amplitude, / = $>/(f>o- At this stage, it is given by a 
Jacobi elliptic cosine, f(x) — cn[x, 1/2]. We have also defined 

k a 1 , 

K = 7m ■ < = T < 69 » 

where k denotes comoving momenta. 



7 In fact, the actual eigenmode equation J68D involves a term — i n the effective frequency. For the background inflaton oscillating around 

the minimum of the quartic potential, the average value of — is zero. However, immediately after inflation, this term gives some residual 
contribution to the effective frequency, which slightly affects the solutions for the eigenmodes and for the background oscillations. 
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During preheating, the modes are amplified by parametric resonance as the inflaton <f> oscillates around the minimum 
of the potential. The key parameter here is the resonance parameter q. Since the coupling A is very small, we may 
expect that q » 1. The modes Xk are then amplified in small intervals of time, when the inflaton crosses zero. 
Parametric resonance amplifies only infrared modes belonging to resonance bands. For large q, the typical momenta 
&;* and widths Afc of the resonance are bounded by [36| 

K, A/c ~ \/A> (70) 

Note that the inflaton fluctuations &k obey an equation similar to the one for Xk, but with q = 3. The resonance 
is not very efficient in this case, so that generally the x-modes are amplified first. However, as soon as the x-field 
reaches high occupation numbers, it excites inhomogeneous inflaton fluctuations through its interaction with cj), again 
with very large amplitude. 

When the amplitude of the fields is sufficiently large, the dynamics becomes fully non-linear, but it may be inves- 
tigated with the help of lattice simulations [33], [38| . This is the stage we will mainly consider in the next section 
when we present our numerical results. Although the details of preheating depend on the particular model considered, 
the non-linear evolution exhibits several generic features [291. At the end of preheating, there is a violent stage of 
rescattering, with a restructuring of the scalar field spectra[37| . [39j]. Around the same time, large bubble- like field 
inhomogcncitics, amplified during preheating, collide to form smaller structures [7j. Finally, there is a much longer 
regime of turbulent interactions, during which the spectra slowly propagate towards thermal equilibrium. 

B. Numerical Results 

We computed the gravitational radiation emitted in the model ([66]) for several values of the coupling constant 
g 2 . The abundance of gravity wave energy density today and its spectrum per logarithmic frequency interval were 
calculated numerically from Eqs. (|47|) and (|28|) . as outlined in section fV Al For the model (|66|) . the universe is radiation 
dominated soon after the end of inflation, so that a 4 pj in Eq. (I47p is very close to the initial energy density at the 
end of inflation (where the scale factor is normalised to unity). The precise value was taken from the simulations: 
a 4 pj ~ 1.15 A0q/4, for all the cases we considered. 

We first present the results of simulations for q = g 2 /X = 120. Fig. [4] shows the gravity wave spectrum (I43[) . 
accumulated up to the time Xf = 240 (t/ = 24O/(-\/A0o) in ([2"5|) ), i.e. well after the end of the preheating stage. 
Each of these spectra were calculated by averaging over different directions in k-space as described in IV Al but the 
directions used were different for the two plots. The two plots were also obtained from simulations with different box 
sizes. Fig. 2] thus confirms that the results are isotropic and independent of the box volume. For this model the peak 
amplitude was h 2 f2 gw ~ 3 x 10~ n with a peak frequency / ~ 7 x 10 7 Hz. The infrared tail of the spectrum at the 
end of the simulation varied approximately as f2 gw oc /. Of course, the spectrum is expected to decay much faster 
at lower frequencies corresponding to modes produced outside the Hubble radius. For this model, the Hubble rate 
dropped from H ~ 2.5 x 10~ 8 Mpi at the end of inflation to H ~ 5 x 10 -11 Mpi at the beginning of preheating, and 
then slowly decreased up to H ~ 10~ 12 Mpi at the end of the simulation. For k/a = H at the time of production, 
this gives a frequency today: / ~ 10 7 Hz, / ~5x 10 5 Hz and / ~ 5 x 10 4 Hz respectively. In particular, a physical 
momentum of the order of the inverse Hubble radius when the main part of the spectrum is produced corresponds 
to a frequency today / ~ 10 5 Hz. All the modes shown in Fig. © were thus well inside the Hubble radius at the 
time of production. This also means that any gravitational waves produced from preheating in this model would be 
at frequencies well above the range of LIGO/VIRGO (to say nothing of LISA or DECIGO/BBO). 

Figure 0] can be directly compared with Figure 1 of [l|| which corresponds to the same parameters and the same 
frequency range. The spectrum obtained in Ref. [l9| has a peak amplitude of h 2 fi gw ~ 10™ 9 at the frequency 
/ ~ 10 8 3 -ffz. In our case, the peak of the spectrum occurs at a lower frequency, and its amplitude is smaller by factor 
of order 1/30. The spectrum we obtain is broader, and in particular it decreases more slowly in the infrared. In the 
next section we will confirm our results analytically. 

We now consider in more detail the role of the different stages of the scalar field evolution on the production of 
gravity waves. Figure [5] shows how the total energy density in gravity waves ([42| is accumulated with time. For 
comparison, we plot on the same figure the evolution with time of the total particle number density of \ arL d <f>. The 
first stage on this plot, up to Xf ~ 90, corresponds to preheating by parametric resonance. During this stage, gravity 
waves are produced in very short intervals where the scalar fields are non-abiatically amplified. A detailed discussion 
of this regime will be given in the next sub-section. The resulting gravity wave energy density grows with twice the 
exponent of the scalar field number density, p gw oc n 2 ot . At the end of preheating, during the short, highly non-linear 
stage of "rescattering", the energy density in gravity waves continues to increase significantly up to i/ ~ 150, when 
the number density becomes approximately constant. For different values of the (/-parameter, this rescattering stage 
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FIG. 4: Spectrum of gravity wave energy density in phys- 
ical variables today (|43[) . accumulated up to the time 
x f = 240, for the model JBBJ with q = 120. The 2 spectra 
were obtained from simulations with different box sizes, 
and averaged over different directions in k-space. 



FIG. 5: The thick curve shows the total energy density 
in gravity waves (|42[) accumulated up to the time Xf , as a 
function of Xf. The thin curve shows the evolution with 
time of the total particles number density, n to t = n x +n,p, 
rescaled to fit on the same figure. 



may be much shorter, and hardly distinguishable form the end of preheating. Finally, the gravity wave energy density 
still slightly increases during the early stage of the subsequent turbulent regime, but at a much lower rate. As shown 
in Figure [5j the exponential rate of gravity wave production is maximal during parametric resonance, but the main 
part of the final energy density in gravity waves is produced during the "bubbly" stage, a short intermediate stage 
between parametric resonance and the onset of the turbulent regime. 




K 



FIG. 6: Spectrum (|71[) of the gravity wave energy density, 
accumulated up to different times function of the 

comoving momentum k (in units of Xcj)o). The spectra are 
shown from Xf = 90 to xf — 240 with spacing Axf = 10. 



FIG. 7: Measure of the (unnormalised) total energy den- 
sity in the two scalar fields per logarithmic momentum 
interval at different moments of time. The same times 
as in Fig. © are shown, the spectra moving towards UV 
from x — 90 to x = 240 with spacing Aa; = 10. 



Let us now discuss how the gravity wave spectrum builds up with time. Here we will consider the spectrum of the 
energy density in gravity waves around the time of their production 



1_ dp gw \ ^ S k {T f ) 
ptot d\nk) p ~ pj 



= ( - ^ ) - ^ 



as a function of the comoving momentum k. Remember that a'j pj ~ 1.15 A0g/4. The spectrum (|71l) . accumulated 
up to different times r/, is shown in Fig. |6]). For a better view, we only plot the spectra starting from the end of 
preheating up to the end of the simulation, because the final spectrum is mainly determined by this stage (see Fig. 1101 
of the next sub-section for a plot of the spectrum during the preheating stage) . For comparison, we show in Fig. [7] 
the spectra of the quantity fc 3 cj£ n£ + fc 3 uif. nf. at the same moments of time, where n k and u)k are the occupation 
number and the frequency of the scalar modes. This quantity gives a rough measure of the total energy density 
in the scalar field fluctuations per logarithmic momentum interval (provided lj^ is changing nearly adiabatically) , 
in analogy with (J7TJ) for the gravitational waves. During preheating, the scalar field spectra are strongly peaked 
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around the resonant momenta inside the resonance bands. The gravity wave spectrum has a peak amplitude growing 
exponentially up to (fig W )p ~ 1CT 9 at xj ~ 90. The subsequent evolution is shown in Figs. The peak of the 
scalar field spectrum is initially in the infrared, around the resonant momentum of preheating, and then quickly moves 
towards the ultraviolet as a result of rescattering. The gravity wave spectrum increases significantly during this stage, 
but in contrast to the scalar fields, the frequency of its peak stays approximately the same. The maximal gravity 
wave amplitude ((f2 gw ) p ~ 5 x 10~ 6 ) is reached around Xf ~ 150. Later on, both spectra slowly propagate towards 
higher frequencies, which slightly increases the total energy density in gravitational waves. On the other hand, the 
peak amplitude and the infrared part of the gravity wave spectrum stay unchanged after Xf ~ 150. (We checked this 
with simulations extending further in the turbulent regime). Their final value may thus be efficiently determined by 
finite-time simulations. 

Now we show the results of simulations for different values of the coupling constant g 2 and hence of the resonance 
parameter q = g 2 /X. They were all performed on lattices with 256 3 points, up to the same final time Xf = 240. 
The spectra of gravity wave energy density accumulated up to different times are shown on the left panels of Fig. [8l 
for q = 1.2, 120, 128 and 1130. On the right panel we show for each of these values of q the evolution with time 
of the spectrum of K 2 \Xk \ + K 2 |$fe| for the two scalar fields. The first observation is that the peak amplitude in 
gravitational waves depends only mildly on the value of q for the cases that we considered. In particular, we don't 
observe the dependence in Q gvJ oc X/g 2 derived in [l8j and quoted in (l9| . The maximal amplitude in gravity waves 
that we observed was for q = 1.2, in which case the spectrum today had a peak amplitude of h 2 f2 gw ~ 5 x 10~ 10 at 
a frequency of order 5 x 10 6 Hz. 

The evolution of the gravity wave spectrum is qualitatively similar to what we observed above. The IR part of the 
spectrum is essentially produced around the end of preheating, when the peak in the scalar fields spectra becomes 
maximal and is quickly shifted towards the UV. Later on, the UV part of tt gw (k) increases, for instance with new 
peaks appearing at high frequencies in the case q = 128 and q = 1130, but the peak amplitude at low frequency 
stays basically unchanged. The frequency of this peak of the gravity wave spectrum depends directly on the resonant 
momentum fc* amplified by preheating (corresponding to the initial peak in the scalar fields' spectra), instead for 
instance of the frequency of the scalar fields' modes (u>k oc *Jq). The corresponding frequency of the gravity waves 
spectrum's peak today is given by l[17|) 

/* ~ fc * 4 x 10 10 Hz ~ K* A 1/4 5 x 10 10 Hz ~ Jf*2xl0 7 Hz (72) 
a i Pj 

where we have used a 4 pj ~ 1.15A0q/4 and fc* = \/A<^o-K*- We took A = 10~ 14 is the last equality. The typical 
momenta amplified by parametric resonance tend to increase as q 1 ^ 4 , see (I70|) . but their precise value is a non- 
monotonic function of q. For instance, for q = 128 (third line from the top in Fig. j8])), we have g 2 /X — 2l 2 with I 
an integer, and in this case fc* is shifted towards zero |36]. Accordingly, the maximal amplitude of the gravity wave 
spectrum in the IR is shifted towards smaller frequencies compared to the case q — 120. For q — 128, the resonance 
is also more efficient. As a consequence, preheating ends earlier and the simulation extends more into the turbulent 
regime. 

In configuration space [7j , the profile of the field inhomogcncities of x and <fr depends directly on the peak momentum 
and width of the spectra k 2 \xk\ an d k 2 \<j>k\- During preheating, the exponential growth of the scalar fields' spectra 
with a sharp peak at the (comoving) resonant momentum fc* corresponds to the amplification of bubble-like fields 
inhomogeneities with a characteristic (physical) size i?* ~ a/k*. At the end of preheating, the broadening of the 
scalar fields' spectra and the shift of their peaks towards the UV corresponds to the fragmentation of these spatial 
fields' fluctutations into smaller structures (see sub-section IIV Bp . Eq. (|53"j) estimates the contribution of this bubbly 
stage to the energy density in gravity waves at the time of production, i.e. shortly after the end of preheating. At 
that time, the total energy density in gravity waves is given by the peak of the spectrum, so this gives 

eu-tX-^W* (73) 

where we took i?* = a/fc* and a p is the value of the scale factor when the peak of the spectrum is reached. We 
factored out the product a 2 H because it is constant at that time and had the same value for the different values of 
q. In all our simulations, the last factor in Eq. (|T3"]) was: (a 2 H) /A^q — 0.28. We checked Eqs. (|5"3"1) . (f?3"j) with our 
numerical results. The dependence on l/a 2 is expected on general grounds because, for the model ([55]) . T? T in l[2"8"]) 
dilutes as l/a 2 with the expansion, giving an overall factor of l/a 2 in Sfe. The peak amplitude of f2 gw is reached when 
the rescaled scalar modes Xk and $>k stabilize to their maximal value. Gravity waves produced after that time are 
further suppressed by the expansion. We also observed the dependence in 1/fc 2 to the extent that we could measure 
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FIG. 8: Evolution with time of the gravity wave energy density spectrum (left panels) and of the scalar field spectrum (right 
panel) for different values of q = g 2 /\. From top to bottom, q = 1.2, 120, 128 and 1130. The left panels show the spectra 
(Q gw)„, see (|71[) . accumulated up to different times Xf, as a function of the comoving momentum K — fc/(\/A</>o). The right 
panels show the total scalar field spectrum 7f 2 (|Xfc| + |$fe|)/0o at different times, as a function of K. 



this with our limited number of simulations. The value of (f2* w )^ tends to decrease with q, for instance from about 

3 x 10~ 5 for q = 1.2 to 3 X 10 -6 for q = 120. However, its precise value varies with q in a non-monotomic way, being 
for instance higher for q = 128 than for q — 120. This is in agreement with the fact that the resonant momentum fc* 
is lower for q = 128. In that case, preheating is also more rapid so that a p is lower. All in all, from our numerical 
results we find that Eqs. (|53|) and (|73|) are satisfied with a ~ 0.15. The precise value may not be significant - because 
we considered only a limited range of values of q and because it is difficult to determine precisely the value of the 
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scale factor "at the time when the peak is reached" - but the order of magnitude is correct. 

To summarize, for the model (1661) . the numerical results for the frequency and the amplitude of the gravity waves 
spectrum are well described by Eqs. (|72l l73|) . They depend very mildly on the details of the rescattering stage, but 
mainly on the characteristic momentum &;» amplified during preheating. The precise value of fc» is a non-monotomic 
function of the parameters, but it may be calculated analytically for any values of the coupling constants. 

C. Analytical Check for Parametric Resonance 

We now study in more detail the production of gravitational radiation from scalar fields experiencing parametric 
resonance with large resonance parameter, q » 1. We consider the linear stage of preheating, when only one field 
has been significantly amplified, so that Eq. (|5T|) is applicable. As we saw in the previous section, the main part of 
the gravity wave spectrum is generated during the subsequent evolution. However, during the linear stage, analytical 
solutions for the modes may be derived, and this will allow us to check our numerical results in details. 

The energy density in gravity waves (|5ip involves time integrals of the form 

1 = f T dr' cos (k t') a(r') Xp (r') X | k _ P | (r') = f X p (x') X {k _ pl (x') (74) 

Jn Jxi V A(p a{x S 

where we have performed the rescalings (|67)l . (f69|) . The basic observation is that, during preheating, gravitational 
waves are produced in a step-like manner, in short intervals of time only, see Fig.[5J These correspond to the moments 
when the inflaton condensate crosses zero and the adiabaticity condition for the scalar modes is violated. The typical 
evolution of the modes with time is shown in Fig. [9] Away from the zeros of the inflaton, the modes oscillate with large 
frequency (oj^ ~ \M)i an d their amplitude evolves adiabatically. In contrast, their amplitude increases sharply when 
the inflaton vanishes (for the modes belonging to the resonance bands) . Only these short intervals of time contribute 
significantly to the time integral in (fT4"|) . Indeed, according to the no-go theorem of sub-section pi Bp . no gravitational 
radiation is emitted by a bath of scalar modes whose frequencies evolve adiabatically with time. The discussion is 
more involved in the case of parametric resonance, because the adiabatic regime then lasts only for a limited period 
of time, between two zeros of the inflaton 8 . But we will confirm below (see Appendix) that it is sufficient to consider 
the vicinity of the moments of time when the inflaton vanishes. 



The situation here depends on the particular channel of interaction considered, namely on the relative signs in the phase of dl9| . as in 
the "particle-like" interpretation of rescattering discussed in [g]. In the analog of l|19|l. the time integral is now performed on a finite 
interval between 2 zeros of the inflaton. Some terms involve a phase of the form 0(x) = J x dx' u> p (x') + J x dx' u>k— — Kx. Because 
u) p ~ y/q, the phase then oscillates very quickly and the time integral reduces to its contribution where 0(x) is stationary. This gives 
the condition of energy conservation for the annihilation of two x-particles into one graviton, which is not possible. On the other hand, 
when the terms in ui p (x') and p (x') in the phase 0(x) have opposite signs, the large terms in <Jq cancel each other and there is no 
direct analogy with energy conservation between interacting particles (the stationary phase approximation is then applicable only for 
large K compared to half the period of the inflaton oscillations). In this case, the interpretation in terms of interacting wave packets is 
more appropriate. 
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FIG. 9: Typical behavior of the modes amplified by parametric resonance. The left figure shows the evolution with time of 
the real part of Xk(x) (see (|67[l ) for k inside a resonance band, here for q = 5050. The right figure shows the evolution of the 
corresponding occupation number, together with the inflaton zero-mode. The 2 figures are taken from j3f|. 



Consider the evolution of the modes around the j zero of the inflaton, at x — Xj. They correspond to waves 
scattered on a parabolic potential @, (36|. Around its zeros (for \x — Xj\ < 1), the inflaton is approximately linear, 
$ oc (x — Xj), and Eq. ([68]) reduces to 



(fXp 
dx 2 



(P 2 + \{x-Xj) 2 ) X p ^0 



(75) 



where P is defined in terms of the comoving momentum p as in (|69[) . The solution is given in terms of parabolic 
cylinder functions W 



X p (x) = A p W 



+ B p W 



P 1 



, -(2 9 ) 1 /4 {x - Xj ) 



(76) 



where A J p and B p are constant coefficients. For x < Xj in the adiabatic regime (for \x 
be solved by the WKB approximation 



X p (x) 



y / 2uj p (x) 



-i J xo 3 Lj p (x') dx' _j_ 



1% 



y / 2uJ p (x) 



a i f x g ui p (x') dx' 



> g" 1/4 ), Eq. ((681) may 



(77) 



where a? p and (3 p are constant (Bogoliubov-like) coefficients. 

In the Appendix we show how to match both regimes ([76]) and (fTTj) . This allows us to follow the whole evolution 
of the modes with time, and therefore to calculate the integral (|74[) involved in the gravity wave energy density. The 
result for the gravity wave spectrum (|27|) , (|5Tj) after the j th zero of the inflaton is then given by 



S: 



(VA0 o ) 6 m } 
a?{x 3 



TT 3 q 



K 



J ^ du (1 - u 2 ) 2 J dP P 6 C P e, k _p| n p +1 nj+^p, (Ic 2 + Is 2 



(78) 



where n p +1 are the occupation numbers for \ i n the adiabatic regime after the j th zero of the inflaton, £ p is defined 

in ([M]) . Ic is defined in (|100p and Is is defined as in (|100[) with the cosine replaced by a sine. Here u = cos(k,p), so 
that |k- p| 2 = k 2 +p 2 - 2kpu. 

We used Eq. ([75| to test our lattice calculation of the gravity wave spectrum. In principle, it is possible to use 
analytical expressions for the occupation numbers given in [361 ] , but it would be difficult to perform the integrals over 
u and P analytically. Instead, we take an interpolation of the occupation numbers output by the lattice simulation, 
and perform the integrals in (|78| numerically. We then compare the result to the lattice calculation in the linear stage 
of preheating, when only the field has been significantly amplified. Fig. [10] shows the resulting spectrum of gravity 
wave energy density. We see that the lattice calculation reproduces Eq. ([78]) very accurately. The spectra resulting 
from the lattice calculation, which correspond to a particular realisation of the initial quantum fluctuations for the 
scalar field, oscillate around the ensemble averaged spectrum obtained from Eq. ([78]) . 
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FIG. 10: Spectrum of energy density in gravity waves at the time of their production, (|71fl . as a function of their comoving 
wave-number k (in units of yX(j)o). The 2 thin lines correspond to the spectrum averaged over different directions in k-space, 
as obtained from the lattice calculation. The thick line corresponds to the spectrum obtained from the analytical formula (1781) . 
The spectra are calculated in the model (166[l with q = 128 during the linear stage of preheating (xf ~ 50 in (|28l) ~). 



In Eq. (|78p , the p-dependance comes essentially from the factor in P 6 n p ?T.|k-p| > which is maximal for p ~ |k — p| ~ 
p*. For k << p*, the two integrals in (|75|) do not depend on k. For K/^flq < few x 2ir, the cosine in (| 1 00[) may go out 
of the time integral, and it cancels with the corresponding sine coming from Is. Therefore, at small wave-numbers, 
the only fc-dependence in Eq. (|78p comes from the factor K 3 in front of the integrals, which determines the infrared 
slope of the spectrum in Fig. [TDJ 



VII. EXPECTATIONS FOR PREHEATING AFTER HYBRID INFLATION 



In this section, we provide preliminary estimates for the gravity waves emitted from preheating after hybrid inflation. 
The mechanism for preheating is different than for the model considered in the previous section, but the resulting 
fragmentation of bubble- like field inhomogeneities is qualitatively similar. We will therefore extrapolate to this case 
the results obtained before for the peak frequency and amplitude of the gravity wave spectrum emitted from the 
bubbly stage. We will analyze this case in greater detail in a subsequent publication (4lj |. 

In the model of chaotic inflation that we considered in the previous section, the energy density during preheating is 
too high for the resulting gravity waves to fall into the frequency range accessible to interferometric direct detection 
experiments (from / ~ 10~ 4 Hz corresponding to the lower limit of LISA, to / ~ 10 3 Hz corresponding to the upper 
limit of LIGO). Indeed, assuming radiation domination at the time of production, the gravity wave frequency today 
(45) may be re-written 



kp Pp 
~H~p 10 8 GeV 



1Hz 



(79) 



where H p and p p are the Hubble rate and the total energy density at the time of gravity waves emission, and k p 
is their physical wave-number at that time. Since we can only expect a significant amplitude in gravity waves for 
wavelengths inside the Hubble radius, k p > H p , it is clear from Eq. (f79")) that models with low energy scales are more 
interesting from an observational perspective. Note however that the peak of the gravity waves spectrum generally 
occurs at higher frequencies, k p >> H p . 

In hybrid models, the inflationary energy scale is typically a free parameter ; it could just as easily be of the order 
of the GUT scale or of the order of electroweak scale. At the end of hybrid inflation, the inflaton decays via a spinodal 
instability of the inhomogcncous modes accompanying the symmetry breaking, in the process of tachyonic preheating 
This also involves inhomogeneous spatial structures, whose typical size i?* depends on the model considered. 
Based on Eqs. (|72|) and (f73|) . we estimate the present-day frequency and amplitude associated with the peak of the 
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gravity wave spectrum emitted from this bubbly stage as 



R*p P ' 



: n* gw ~ io- 6 (iz, # p ) 2 (8i) 



Hybrid inflation should generally end quickly compared to the expansion rate of the universe, the so-called waterfall 
constraint. It follows that for practical purposes, one may usually neglect the expansion of the universe during 
preheating. The energy density p p above is then just given by the energy scale during inflation, p p — Vi n f, and H p is 
the corresponding Hubble rate. Note that the waterfall constraint implies typically i?* > 1/H p , so this already places 
some lower limit on /* and some upper limit on f2* w . 

We will estimate and (f8"Tj) as a function of the parameters for two particular models of tachyonic preheating. 
A prototype of hybrid inflation model is given by the potential 



V=\(* 2 -vi) + ^4> 2 v 2 + ^<f (82) 

For <f> > (f) c — ^ v, the fields have positive effective mass squared, and the potential has a valley at a = 0. Inflation 
occurs while (f> decreases slowly in this valley. The energy scale during inflation is usually dominated by the false 
vacuum contribution, V\ n f = j v 4 . When <j) reaches the bifurcation point <j> c , inflation ends and the fields roll rapidly 
towards the global minimum at <f> — and \a\ = v. To avoid the production of domain walls, one usually considers a 
complex field a. 

When reaches <j> c , the curvature of the effective potential becomes negative and field inhomogeneities are amplified 
by the tachyonic instability. The details of this process, in particular the typical momenta fc* that are amplified, 
depend on the particular trajectory followed by the fields <fi an d which in turn depends on the values of the 
coupling constants and on the initial inflaton velocity at (f> = <p c . Here we will only consider two particular trajectories 
in field space, which lead to different predictions for the gravity waves spectrum. Preheating in these two cases has 
been studied in detail in [f|, @, and we will rely on their results for our estimations. 

For g 2 >> A, or sufficient inflaton velocity at <f> = 4> c , the inflaton overshoots the bifurcation point and the field a 
falls down only when when << <j> c . In this case, the effective potential for a is dominated by [8J] 

V=\*-\#f + \d* (83) 



The typical momenta amplified by preheating in this model vary as V~\v/ y In (^p) ■ From the numerial results 
presented in [9(, we estimate the typical size of the inhomogeneous structures of a when they fragment as 



The resulting frequency (f8T)|) and amplitude (|5Tj) associated with the peak of the gravity wave spectrum are then given 
by 

/»~ ai/ ; iqiohz (85) 



and 



^~2*10-Mn(^)(-iL) (86) 

Note that /» depends on the energy scale only through the coupling constant A, and that this dependence is very weak, 
essentially varying as A 1 / 4 . The peak of the gravity wave spectrum may fall into the range accessible to interferometric 
direct detection experiments, /* < 10 2 Hz, only for extremely small values of the coupling constant, A < 10~ 28 . 

Another special case corresponds to g 2 = 2 A in (152")) and sufficiently low inflaton velocity at the critical point 
(f> = (f> c . In this case, the effective potential has the same structure as the effective potential of SUSY-inspired F-term 
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inflation, and the fields <f> and cr fall down along a simple linear trajectory, see e.g. [4(j. Along this trajectory, the 
effective potential for the field a takes the form (up to the redefinition A — > A/3) 

Preheating in this model has been investigated in detail in [9| . The typical momenta amplified vary as fc* ~ . The 
process is initially dominated by fluctuations with momenta somewhat greater than k*, whose initial amplitude is 
more suppressed, but which grow faster. It follows that the bubble-like field inhomogeneities initially have a typical 
size i?* somewhat smaller than but they are also quite distant from each other and they expand for some time 

before colliding. We estimate the typical size of these spatial structures when they collide from the numerical results 
of 

R* ~ ^- (88) 
Xv 

The resulting frequency (|80[) and amplitude (|8ip associated with the peak of the gravity wave spectrum are then given 
by 

/» ~ A 3/4 10 10 Hz (89) 

and 

Note that this case generally leads to higher amplitude and smaller frequency in gravity waves than the previous 
model (|83p. Requiring the frequency of the peak of the gravity wave spectrum to fall into the range accessible to 
interferometric direct detection experiments, /* < 10 2 Hz, gives A < 2 x 10 This is still a very small value for the 
coupling constant, but it is considerably higher than for the model (|83l) . 

For example, according to Eqs. (|59"|) and (|9H)) . a model with g 2 = 2 A = 2 x 10~ n and v = 10 12 GeV (so that the 
waterfall constraint is marginally satisfied) would lead to /* ~ 50 Hz and h 2 fi* w ~ 2 x 10~ 9 , which could already be 
detectable by LIGO. Note that we only considered here the peak of the gravity wave spectrum, but the infrared tail 
of the spectrum could also be observationally relevant. 

The preliminary estimates made in this section are very simplified and should be supplemented by numerical 
simulations [4l| . We considered only the peak of the gravity wave spectrum emitted during the bubbly stage for two 
particular models, neglecting the amplification of the field <$> or any degree of freedom other than a. In particular, other 
models should be considered, and a significant amount of gravity waves may still be produced during the transition 
from the bubbly stage to the regime of well-developed turbulence 9 . However, we see that even in simple models, the 
results for the gravity wave spectrum may cover a wide range of values. 



VIII. DISCUSSION AND SUMMARY 



In this paper we constructed a general formalism to calculate the spectrum of classical gravitational waves emitted 
from random media of dynamical scalar fields in an expanding universe. In principle this can be applied to any 
cosmological situation where scalar field dynamics is important and potentially interesting for the production of 
gravitational radiation. This includes, for instance, early universe phase transitions at the electroweak energy scale 
or above [l2l |. the decay of the inflaton field in the process of (p)reheating after inflation, and other non-equilibrium 
processes. 

In the paper we advanced the theoretical framework for preheating after inflation, which is accompanied by the 
violent decay of the inflaton and the copious production of bosonic fields coupled to it. Based on our formalism, we 
have developed a numerical code to calculate gravity wave production. 

The formalism and numerics for gravity wave emission in this paper were implemented for models of chaotic inflation, 
and in particular for chaotic inflation with a quartic potential. This case served as a convenient ground to contrast 



In the chaotic inflation model studied in the previous section, the energy density in gravity waves produced during this transition stage 
is suppressed by the expansion of the universe, but we don't expect such a suppression for hybrid inflation models. 
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our numerical results with previous ones, and to test them analytically in great detail. The dominant contribution to 
gravity wave emission comes from the non-linear "bubbly" stage at the end of preheating. Our numerical results for 
the present-day frequency and amplitude of the resulting spectrum's peak are well described by the simple formulas 

4 x 1 10 Hz 

f* and ^w~10~ 6 (R*H p f (91) 

R* Pp 

where H p and p p are the Hubble rate and total energy density during the bubbly stage. The characteristic (physical) 
size i?* ~ a/k* of the bubble- like field inhomogeneities depends on the typical (comoving) momentum fc* amplified by 
preheating. Up to factors depending on the expansion of the universe, /* and h 2 f2* w are then known functions of the 
parameters for different models. The present day frequency /» depends on the particular chaotic model considered, 
but it remains too high, above 10 6 Hz, to be detected by currently planned experiments. The fraction of energy density 
in gravity waves also varies from one chaotic model to another. It is generally lower for the model m 2 4> 2 + g 2 (f> 2 x 2 
than for the model with A <^ 4 inflation, because the typical momenta amplified fc* ~ y/gm<&o are higher in that case. 
The results for p gw /ptot may also be different in models with trilinear or non-renormalisable interactions (3lj |. where 
the dynamics ocurs faster so that (i?* H) is less diluted by the expansion of the universe. 

Based on Eq. (|9ip . we also estimated the frequency and the amplitude of the peak of the gravity- wave spectrum 
produced from the bubbly stage in two particular models of preheating after hybrid inflation. Because there are more 
free parameters in this case, the results may cover a much wider range of values. We found that in principle, for models 
with very small coupling constants, the resulting gravity waves spectrum may be relevant already for LIGO/VIRGO. 
Gravity waves from preheating after hybrid inflation certainly deserve further investigations. 

Besides gravity wave interferometers probing typical frequencies which range from 10 Hz (LISA) up to 10 3 Hz 
(LIGO), there are several bar or spherical resonant detectors (see e.g. (42[) operating in the kHz range. Experiments 
have also been proposed at higher frequencies, up to 100 Mhz 43]. However, the sensitivity to h 2 il gw drops when 
the frequency increases, and the gravity waves produced from preheating after chaotic inflation lie outside the range 
currently accessible to these experiments. 

The frequencies of gravitational waves emitted from preheating could be naturally redshifted to lower values if 
there was some intermediate matter dominated stage between preheating and the radiation dominated era. If the 
scale factor expands by a factor a m during this stage, then the frequencies of gravity waves emitted before that stage 
are decreased by a factor of 1/a™ 4 . However, the fraction of energy density in gravitational radiation is then diluted 
by a factor of l/a m . 

We only considered the dynamics of two coupled fields during preheating, the inflaton and another scalar field. 
Clearly, the production of gravity waves will increase with the number of degrees of freedom. In particular, non-scalar 
bosonic degrees of freedom excited during preheating will change the results for gravity wave production during the 
turbulent evolution towards thermal equilibrium. Indeed, we found that neither established Kolmogorov turbulence 
for scalars nor scalars in thermal equilibrium contribute to p gvJ . The situation will be different for gauge bosons and 
photons. It will be interesting to address the issue of gravity waves emission from the vector fields excited during 
preheating. 

IX. ACKNOWLEDGEMENTS 

We would like to thank Latham Boyle and Andrei Frolov for stimulating discussions, and Richard Easther, Daniel 
Figueroa, Juan Garcia-Bellido and Eugene Lim for useful discussions and correspondences during different stages of 
the project. A.B. and G.F. were supported by NSF grant PHY-0456631. L.K. was supported by NSERC and CIFAR. 

X. APPENDIX 

In this Appendix, we derive the evolution with time of the modes X p of sub-section lVI Cl during parametric resonance 
in the model (|66p . This allows us to calculate analytically the gravity wave energy density spectrum emitted during 
this stage. 

Around the zeros of the inflaton field, the modes X p are given by Eq. ([76]) . In between the zeros of the inflaton, 
they are given by the WKB solution (|76| . In this regime, the occupation numbers are constant and well defined. 
With our definition of (3 J above, they are given by 



n p = V\<j> \ti\ 2 



(92) 
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where the prefactor comes from the fact that u>k is the rescaled frequency, see (|68|) and (|69|) . Note that X p and /K 
have dimension mass -1 / 2 , as the modes x P in (|49|) . so that n^. is dimensionless. 

Both solutions (|76f and (|77| overlap on an interval of cc of order q -1 / 4 < |x — Xj| < 1. Matching the four coefficients 
gives 



B 3 = 
p 



(2?) 



-1/8 



j4 e -i(T/4+0j-C P /2) + el (7r/4+0J-Cp/2) 



4 = £ Wf (2q)' 1 / 8 a p e -^l^l-^IV + % e *W4+«J-C,/2) 



(93) 



where 



C P = fl + e - 2 - p2 /v^ 



1/2 



3 -7rP 2 /v^9 



(94) 



The phases 8 J p and £ p are given in [y], [36j |. Their precise value will not be important for our purpose. For x > Xj, the 
modes are again given by the adiabatic solution (I77|) with new coefficients a p +1 and (3 3 p +1 . Matching this new solution 
with (|76p gives the evolution of a p , (3 P and n p around the j zero of the inflaton. For large occupation numbers, n J p +1 
at x > Xj is related to n J p at x < Xj by @ 



2 e - 2 ^ 2 /v^ _ 2 S in 6£ e - - p2 /v^ 



(95) 



where 6£ = 2 6^ — ( p + arg/3^ — arga p . 

Because the expansion of the universe may be scaled out of the equations for the model (|66p , the resonance involves 
separate stability and instability bands for the momenta p, which depend on the phase Q J p . Only the modes belonging 
to the instability bands are amplified and contribute to the p- integral in the gravity wave energy density (|51|) . The 



maximal occupation number is reached for sinO^ = —1 in (|95p . corresponding to constructive interference, and it is 
a good approximation to take this value for all the modes under consideration. In this case, Eq. (|93|) simplifies to 



,i(v/4+8$-Cp/2)) 



and A> = 



(96) 



where we have used the fact that a p and (3 p differ only by a phase in the limit of large occupation numbers. Finally, 
using (|92|) and (f95|) with sinG^, = —1, we may relate B 3 p to the occupation numbers before and after x = Xj 



\Bl\ 2 = (2?) 



-1/4 



2£ P 



£ P VA0 o (2<7) 1/4 VXcf>o 



(97) 



We are now ready to calculate the time integral (|74[) . Since the amplitude of the modes is increasing exponentially 
with time, it is sufficient to perform the integral over the latest inflaton oscillation only. The solution (|76|) has a peak 



amplitude at x = Xj, just after x = Xj. It is accurate over an interval of x of order Ax = few x q^ 1 ^ 4 , around ; 
and differs from the exact solution only far inside the adiabatic regime. The amplitude of the modes at that time is 
lower than their amplitude at x = Xj by a factor of order g -1 / 8 , as may be seen from Eqs. (|76|) . (|77jl and (|93|) . This 
gives a relative contribution of order q~ x / 2 to the gravity wave energy density (|5ip . which is negligible for q » 1. 
The integral ([74]) after the j th zero of the inflaton is therefore well approximated by 



P 



+ i 



VA0 o o( 



dx' cos[Kx'] W 



p2 



2q 



-{2q)V\x-Xj) 



W 



IK -PI 



2q 



, -{2q) 1 ' 4 (x - Xj ) 



(98) 



where e = few x q 1 / 4 . The scale factor a(x) goes out of the integral since it is approximately constant on the interval 
[xj — e , Xj — e]. Introducing the variable v — (2q) 1 ^ 4 (x — Xj) and using (|97p . we have 



1 1 " (Va^o) 4 / pC|k - p|np n i k -pi lc 



(99) 



where 



Ic 



Vf 



dv cos 



, — K Xj 



W 



P 2 



2q 



W 



|K-P| 2 



(100) 
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The result does not depend on the precise values of Vi and Vf provided the integration range [vi, Vf] extends sufficiently 
into the adiabatic regime. The functions W depend only slightly on their first argument, which is smaller than one for 
modes inside the instability bands (see ([70])). They are oscillatory for all v, with a peak amplitude at v = v ~ 1.5. The 
amplitude of the oscillations decreases as v moves away from v, and the oscillations at v > v have larger amplitude 
than the ones at v < v. The main contribution to the integral comes from the oscillation centered around v, whose 
width is about Av ~ 4. For a precise calculation, we may also take into account the first oscillation before v and the 
first two oscillations after it, corresponding to Vi ~ —4 and Vf ~ 7. 

Finally, inserting into (|5ip and performing the integral over the azimuthal angle of p, we derive Eq. ([75|) for 
the gravity wave energy density spectrum after the j th zero of the inflaton. 
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